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We consider diffusion of vibrations in 3d random harmonic lattices with translational invariance. 
Above some frequency u>ik, corresponding to the loffe-Regel crossover, notion of phonons becomes 
ill defined. They cannot propagate through the lattice and transfer energy. Nevertheless most of 
the vibrations in this range are not localized. We show that they are similar to diffusons introduced 
by Allen, Feldman et al., Phil. Mag. B 79, 1715 (1999) to describe heat transport in glasses. The 
crossover frequency ojir is close to the position of the boson peak. Changing strength of disorder 
we can vary from zero value (when rigidity is zero and there are no phonons in the lattice) up 
to a typical frequency in the system. Above cjir the energy in the lattice is transferred by means of 
diffusion of vibrational excitations. We calculated the diffusivity of the modes D{ijj) using both the 
direct numerical solution of Newton equations and the formula of Edwards and Thouless. It is nearly 
a constant above cjir and goes to zero at the localization threshold. We show that apart from the 
diffusion of energy, the diffusion of particle displacements in the lattice takes place as well. Above 
wiR a displacement structure factor ^(q, to) coincides well with a structure factor of random walk 
on the lattice. As a result the vibrational line width r{q) = Duq^ where Du is a diffusion coefficient 
of particle displacements. Our findings may have important consequence for the interpretation of 
experimental data on inelastic x-ray scattering and mechanisms of heat transfer in glasses. 



PACS numbers: 63.50.-x,65.60.-l-a,78.70.Ck 
I. INTRODUCTION 

Propagation of vibrational excitations in disordered 
systems is one of the advanced problems in condensed 
matter physics. In particular, transport mediated by 
these excitations is responsible for the thermal conductiv- 
ity of amorphous dielectrics (glasses). However mecha- 
nisms of heat transfer in glasses above the plateau region 
are still poorly understood. 

At low temperatures below 1 K the low frequency plane 
long wave acoustical phonons are well defined excitations 
which transfer the heat in glasses. At these temperatures 
the thermal conductivity x{T) oc and is controlled by 
a resonant scattering of phonons on two-level systems 
(TLSP^l. Between 4K and 20 K the thermal conductiv- 
ity >f(r) saturates and displays a well known platearP. 
As was shown irP it can be explained by resonant scat- 
tering of phonons by quasilocal vibrations (QLV). The 
QLV, together with TLS and phonons are vibrational 
excitations responsible for many universal properties of 
glassed. 

However, above approximately 20 K the thermal con- 
ductivity rises again and finally saturates on the level 
of one order of magnitude higher, i.e. at temperatures 
about several hundreds KelvirP. As generally believed, 
the origin of this second rise of the thermal conduc- 
tivity (above the plateau) is not related to phonons. 
It was established long agcP"^, that in this tempera- 
ture (frequency) range the mean free path of phonons 
I becomes of the order of their wave length A (or even 
smaller, of the order of interatomic distance). Corre- 



spondingly, the loffe-Regel criterion for phononJiS! he- 
comes violated. The existence of such crossover was con- 
firmed by molecu lar dy namics calculations fo r som e real 
and model glasse^^^ and disordered lattice^^2lll_ 

In the regime of such strong scattering a standard 
concept of plane waves (phonons) with a well defined 
wave vector q becomes inapplicable. The question then 
arises: what physical mechanism is responsible for the 
heat transfer in glasses in this temperature range? The 
numerical simulations show that majority of the vibra- 
tional modes in the corresponding frequency range are 
not locahzecP^Jti^. 

As was shown lower limit of the thermal con- 

ductivity of amorphous solids above 30 K can be cor- 
rectly estimated within the framework of the Einstein's 
modePIl. It was assumed that the mechanism of heat 
transport above the plateau is a random walk of thermal 
energy between clusters of neighboring atoms vibrating 
with random phases. In fact, a diffusion mechanism for 
the heat transfer in this temperature range was proposed. 

At the same time, delocalized vibrations in glasses of a 
new type, different fronr phonons, were introduced. They 
were called diffuson^^^^. These are vibrations spread- 
ing through the system not ballistically, as phonons (on 
distances of the order of mean free path) but by means of 
diffusion. It is an important class of excitations which oc- 
cupy in glasses the dominant part of the spectrum^^'. In 
these papers it was put forward the hypothesis that the 
boundary between phonons and diffusons is determined 
by the loffe-Regel criterion for phonons. Since diffusons 
are delocalized excitations, they may be responsible for 
the thermal conductivity of glasses above the plateau. 
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The similar conclusion was made by the authors o P^ l ^^ l 
They considered the case of strong scattering of phonons 
in disordered lattices with a significant fraction of ran- 
domly located missing sites, but which is still far from 
the percolation threshold. It was shown that, in contrast 
to the electronic case, the lofFe-Regel criterion is inaccu- 
rate in the prediction of phonon localization. Instead of 
localization, the vibrational transport above the loffe- 
Regel threshold becomes diffusive with approximately 
constant energy diffusivity D{ui). The diffusivity was cal- 
culated by numerical solution of the Newton equations 
for particle displacements. Similar calculations but for 
real glasses were done in the papers-'^^ ''^ using molecular 
dynamics methods. 

The diffusons above the loffe-Regel crossover were 
identified also in granular jamm ed sy stems with repul- 
sive forces between the particle^^^'^. They also have 
diffusivity which is independent of frequency uj. It was 
calculated making use of the Kubo-Greenwood formula 
for the thermal conductivity derived in^^. In jammed 
systems the loffe-Regel crossover frequency wjr can vary. 
It is shifted to zero when the system approaches the jam- 
ming transition point. 

Therefore, as we believe, it is important to study prop- 
erties of diffusons systematically in systems where they 
exist. They bring a new physics to our understanding 
of vibrational properties in strongly disordered systems 
and energy /heat transfer in glasses. As is well known 
the number of diffusing physical quantities coincides with 
the number of integrals of motion. Besides energy, we 
have two quantities which are conserved in a disordered 
closed mechanical system which is not affected by ex- 
ternal forces. These are the total momentum and the 
position of the center of inertia. Since they are con- 
nected with each other, it should exist at least one new 
additional diffusion coefficient, different from the energy 
diffusivity D{uj) investigated iiPSttgej show in this 

paper that this new coefficient is a diffusivity of particle 
displacements, D„. 

To study these important properties we should have a 
model, sufficiently simple, to describe all of them. Re- 
cently, using a random matrix approach, we developed a 
simple scalar model of a Zd disordered harmonic lattice 
with translational invariance which, as we believe, can 
represe nt typ ical properties of vibrations in disordered 
systemfPEir jt was shown that for a certain type of 
disorder the rigidity of the system is zero and the usual 
phonons (plane waves) cannot propagate through the lat- 
tice. However almost all vibrations in the lattice are delo- 
calized and therefore can participate in the heat transfer. 
As we will show in the present paper they are similar to 
diffusons introduced irpSl. 

One can continuously bring rigidity and phonons back 
to the lattice increasing the width of the so-called phonon 
It is a region where density of vibrational states 
g{oj) oc cj^. We will show that in this gap phonons are 
well defined plane wave excitations. However above the 
gap they cease to exist, since their mean free path be- 



comes of the order of wave length what corresponds to 
the loffe-Regel crossover. 

Nevertheless vibrations outside the gap remain delocal- 
ized and can participate in the heat transfer. Similar to 
the case of zero rigidity the character of the energy trans- 
fer in this frequency range is diffusive and is drastically 
different from the ballistic phonon mechanism. We be- 
lieve that in this case vibrational properties of our system 
may be similar to the properties of real glasses. There- 
fore we are going to consider in more details properties 
of phonons and diffusons in these disordered lattices. To 
start with, for the sake of clarity of further consideration, 
we outline below the main properties of our model. 

II. A RANDOM MATRIX APPROACH 

The vibrational properties of a mechanical system of 
N particles are determined by the dynamical matrix 
Mij = ^ijl ^Jfn^roj, where $jj is the force constant ma- 
trix and vrii are the particle masses. The matrices M 
and $ are real, symmetric and positive definite matri- 
ces N X N (we consider for simplicity a scalar model). 
The last condition is important. It ensures mechanical 
stability of the system. 

One can always present every real, symme tric an d pos- 
itive definite matrix M in the following fornPHMl 

Af = AA^, or Afy = ^A.feAjfc. (1) 

k 

Here A is some real matrix of a general form (not neces- 
sarily symmetric). And, vice versa, for every real matrix 
A the product AA'^ is always a positive definite symmet- 
ric matrix. 

For a free mechanical system it is necessary to satisfy 
also translation invariance conditions 

^M,,=^M„=0 (2) 

(for simplicity we consider below all masses = 1). It 
ensures that the potential energy of the system depends 
only on the differences of particle displacements. 

We believe that some important properties of glasses 
can be reproduced if we take matrix A as a random one. 
As irP2l, we consider the case of a simple cubic lattice with 
N particles and lattice constant qq = 1. Each particle 
has its unique integer index i which takes values from 1 
to TV. We construct the matrix A as follows. The non- 
diagonal elements Aij (for i ^ j) we take as independent 
random numbers from Gaussian distribution with zero 
mean {Aij) = and unit variance (^fj) = 1 if i-th and 
j-th particles are nearest neighbors. For each particle 
in a simple cubic lattice there are six nearest neighbors. 
Non-diagonal elements Aij and Aji are statistically in- 
dependent from each other (matrix A is non-symmetric) . 
All other non-diagonal elements (for non-nearest neigh- 
bors) Aij = 0. To ensure the translational invariance the 
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diagonal elements An are calculated as follows 



(3) 



Then, according to Eq. ([T]), the Eq. ([2]) will be also met. 
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FIG. 1: The normalized DOS g(i^) for random matrix M = 
AJ^ built on a simple cubic lattice with A*' = 20 x 20 x 
20 particles and averaged over 1000 realizations. Inset: The 
participation ratio P(bj) for N = 10^ (a) and TV = 27^ (b) for 
one realization. 

Fig. [T] shows the normalized density of vibrational 
states (DOS) 5(1^) of matrix M — AA^ in this cubic 
lattice. The periodic boundary conditions were used. As 
follows from the figure, g{uj) is nonzero at a; — )■ 0. In spite 
of the presence of translational invariance we do not see 
the expected phonon modes with their DOS (7ph(w) oc 
for a; — > 0. It means that phonons as plane wave excita- 
tions cannot propagate through the lattice. 

As was shown in'^'* it is because the affine assumptions 
are violated and the macroscopic elasticity theory be- 
comes inapplicable in this case. The average value of the 
Young modulus of the lattice E cx Therefore in the 

thermodynamic limit {N — ^ 00) i? — >■ 0. As a result the 
rigidity of the lattice and sound velocity are also zero. 
This unusual behavior is due to presence of high concen- 
tration of negative springs in the lattice which makes it 
extremely soft. 

To determine whether vibrational modes are localized 
or delocalized we have calculated the participation ratio 



N 



(4) 



Here ei{uj) is i-th particle normalized eigenvector with 
frequency lu. As one can see from the inset of Fig. [T] all 
modes with exception of small high frequency part are 
delocalized. They have f (w) « 0.2 which is independent 
of the system size. This value is close to the theor etical 
value 1/3 for Porter-Thomas distribution of ef (wp^*^. 
We have verified also that the level spacing distribution 
obeys the Wigner-Dyson statisticP^l. it also indicates the 
mode delocalization. As we will show in Section IIVI all 




FIG. 2: The spacial eigenmode structure of random matrix 
M — AA''" for the lowest frequency tJmin in two dimensional 
square lattice 400 x 400. 



these delocalized vibrational modes can be identified as 
diffusons. They spread in the lattice by means of diffu- 
sion. 

To elucidate a spacial structure of the eigenmodes for 
matrix M — AA^ we considered as an example a two 
dimensional square lattice with N = 400 x 400 particles 
and calculated eigenvector ei(a;i„in) (« = l,2,...,iV) for 
the lowest frequency Wmin in the system. The result is 
shown on Fig. [2j Particles with positive and negative 
displacements are shown by white and black dots corre- 
spondingly. As one can see from the figure, the mode is 
delocalized. Its spatial structure is random (fractal) and 
has nothing to do with a plane wave. Similar picture 
takes place in a 3d case. 



III. PHONONS 

To introduce phonons into the picture we should have 
finite rigidity of the lattice. Since a sum of positive def- 
inite matrices is a positive definite matrix then a possi- 
bility is to add to the random matrix AA^ a crystalline 
par " 



M = AA^ + fiMo 



(5) 



Here matrix A is the same random matrix built on a 
3c? simple cubic lattice with translational invariance and 
flo = 1 as in the previous Section. Matrix Mq is a posi- 
tive definite crystal dynamical matrix for the same lattice 
with unit masses, and all spring constants (between the 
nearest neighbors) equal to unity. The tune parameter 
/Lt ^ controls the rigidity of the lattice. 

To find the rigidity (as a function of fj.) we calculated 
numerically the Young modulus E of the lattice with dy- 
namical matrix given by Eq. ([5| for 7^ 0. For that we 
fixed particles on the left hand side of our cubic sample 
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FIG. 3: Young modulus £ as a function of jj, for dynamical 
matrix M — AA'" + /iMo built on a cubic lattice with A'^ = 
100 X 100 X 100 particles (one realization). Black dots are 
calculated values, the line is the best least-square fit. 

with N — particles and displaced all particles on the 
opposite (right hand) side by the unit distance. In other 
two directions we used the periodic boundary conditions. 
Then, solving the system of linear Newton equations, we 
found the new equilibrium positions of all other particles 
in the sample and calculated restoring forces acting on 
the displaced particles on the right boundary. Due to 
randomness of the elastic bonds, the restoring forces are 
also random. Let / be the average restoring force. Then 
the Young modulus E can be calculated as follows 

E-jfL.. (6) 

To avoid confusion, we remind that we are using here a 
scalar version of the elasticity theory. Therefore all forces 
in the lattice are parallel (or antiparallel) to the particle 
displacements. 

The results of these calculations are shown on Fig.[3]for 
cubic sample with N = 10^ particles. As we can see from 
the fit, the Young modulus has a following dependence 
on /i: 

E ^ fi, > 1, (7) 

E = l.5y/JI, (8) 

As a result, for /x 3> 1 we have a usual crystal, where 
disorder is relatively small and relation ([t]) is obvious. 
For <C 1 the disorder is strong. The fluctuations of 
the nondiagonal m atrix elements Mi^j are bigger than 
the averaged valueJ^SIH. In this case Young modulus 
E cx The origin of this behavior is unclear and it 

should be elucidated in future work. But we will support 
below our numerical findings by calculation of the sound 
velocity and the phonon density of states (for small uj) 
and comparison of the latter with total DOS calculated 
numerically for matrix ([5|. 

To calculate the phonon contribution to the DOS at 
small Lo, we need to know the sound velocity v at zero 
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FIG. 4: The normalized DOS g{ijj) for dynamical matrix 
M = AA'^ + fiMo and diflterent /i (0, 0.001, 0.01, 0.1, 1) 
calculated with precise numerical KPM solution for cubic lat- 
tice with A'^ — 200^ (full lines). Straight lines correspond 
to Eq. ( |l0| with sound velocity v = Filled and open 

diamonds correspond to phonon contribution to the DOS be- 
low and above the loffe-Regel crossover frequency ujm corre- 
spondingly (see further text for details). Inset: dependence 

frequency. It is related to the Young modulus in a usual 
way: 

v = Ve (9) 

(since all particle masses = 1 and lattice constant 
ao = 1). Then for the phonon DOS (in the scalar model) 
we have 

5phH = ^^. (10) 

The total DOS 5(w), normalized to unity and cal- 
culated numerically by the kernel polynomial method 
(KPM)'^'* for dynamical matrix ^ for different val- 
ues of /i, is shown on Fig. [4] We see from the figure that 
for finite values of /i 7^ the DOS at low enough frequen- 
cies is proportional to which corresponds to acoustical 
phonon excitations. Thus, introducing finite values of /i, 
we open in the vibrational spectrum a phonon gap. Just 
above this gap the DOS has a sharp maximum at fre- 
quency Wmax which wc will identify with the width of the 
gap. As follows from the figure, the maximum frequency 
increases as Wmax y/JI. 

Since the DOS 5(1^) is normalized to unity, we con- 
clude from the Fig. [4] that vibrations corresponding to 
the maximum for /i 7^ were pushed out from the region 
of small frequencies uj < Wmax for /i 0. As we will 
show further (see Table I), the frequency Wmax is corre- 
lated with position of the boson peak cjf, (the maximum 
in the reduced DOS g{uj)/ijj'^). Therefore appearance of 
the boson peak in disordered systems is not necessarily 
related to the acoustic van Hov e singularity in crystals 
as was proposed recentljD^'^SIlI], 
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The straight lines on the Fig. [4] correspond to the 
phonon DOS (7ph(<^) determined by Eq. (10 1 with the 
sound velocity v = \fE and E calculated from Fig. [s] 
One can see a good agreement of the total g{io) at low fre- 
quencies with the phonon contribution (7ph(w). From that 
we can conclude that at least the low frequency excita- 
tions in the phonon gap are the usual long-wave phonons. 
However actually, as we will show further, nearly all ex- 
citations in the gap up to the frequencies close to w^ax 
belong to phonons but with a nonlinear dispersion law. 



field as follows 
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FIG. 5: Participation ratio for different as a function of to 
for A*' = 27'^ (one realization). The arrows indicate positions 
of a;niax in g(yi) for corresponding values of ^ (see Fig. [4|. 

This conclusion is supported by calculations of the par- 
ticipation ratio -P(i^). It is shown in Fig. [S] for various 
values of /i. For /i 7^ one can clearly distinguish in the 
function P(w) two different frequency regions. As follows 
from Fig. |4j the low frequency part (below Wmax) corre- 
sponds to the phonons. In this range the participation 
ratio increases with decreasing frequency. It is related 
to increase of the phonon mean free path 1(lS) as o; — ?> 
(see Fig. |9]). In the high frequency part (above Wmax) 
P(w) is approximately independent of the frequency and 
coincides with participation ratio for /x = 0. As we will 



show in Section IV it corresponds to diffusons. A similar 
behavior of the participation ratio was found recently in 
2d Lennard- Jones glasses*^. 

To find the phonon dispersion curve (dependence of 
the phonon frequency a; on the wave vector q) and 
phonon mean free path Z(w) we should calculate space 
and time Fourier transform of the particle displacement 
field u{r,t). For that we ascribed to all the particles at 
the initial moment t = random displacements (with 
Gaussian distribution with zero mean and unit variance) 
and zero velocities. Then, numerically solving Newton 
equations (using Runge-Kutta-4 method with the time 
step At = 0.01), we analyzed the particle dynamics at 
t ^ 0. 

Let u{ri, t) be the i-ih. particle displacement as a func- 
tion of particle coordinate and time t. We define the 
displacement structure factor (DSF) of the displacement 



5(q,u;) 



2 



N 



(11) 



For better frequency resolution the upper time limit T 
was taken sufficiently large (T = 3000), while the inte- 
gration time step was chosen as At — 0.01. Since vectors 
Ti in a cubic lattice are discrete, the wave vectors q = q„ 
are also discrete and are defined on the corresponding re- 
ciprocal lattice. For example, for cubic sample Lx Lx L 
and q || (100) direction we have g„ = 2TTn/L where inte- 
ger numbers n are — L/2 ^ n ^ L/2. 

One can show that definition ( |Tl| ) is equivalent to the 
usual expression 



N N 

S{q,uj) =7r^|^e,(a;j) 



Siuj-ujj) (12) 



where it is written in terms of eigenvectors ei{ujj) and 
eigenfrequencies Uj of matrix M. The density of states 
is related to the structure factor by the sum rule 



5(a.) = -^5(q,c.). 



(13) 



According to definition (11) S'(0,w) = since the posi- 
tion of center of inertia is conserved and uiji, t) = 0. 
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FIG. 6: The Lorentz dispersion curves for different wave vec- 
tors q II (100) direction and ^ — 0.1. Closed diamonds cor- 
respond to the calculated values of S(q, lu) and lines are fit- 
ting curves according to Eq. (14 1. The number of particles 
TV — 50^ (one realization). Insets: the Lorentz dispersion 
curves for q = 0.5 and q — 0.75. 

To analyze phonon excitations, we have found the max- 
imum of 5'(q, io) as a function of w for each discrete value 
of q„, for several values of /z. As an example, the results 
for /i = 0.1 and one q direction are shown on Fig. [6] For 
the fitting curves we used the Lorentz distribution 
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{UJ — CJq)^ + (Aw)^ 



(14) 
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FIG. 7: The dependence ojq on q for q || (100) direction for various /i (1, 0.1, 0.01) in a cubic sample with A*' — 50^ (one 
realization). Filled and open diamonds are the maximums of S{q,ij) as a function of u for each discrete value of for 
frequencies below and above the loffe-Regel crossover correspondingly (see text below for details). Solid lines correspond to 
halves of the maximums. Dashed lines show lj = vq linear dependence with sound velocity v — V~E. Horizontal dotted lines 
correspond to the maximum frequency ojmax in g(i^) (taken from Fig. [4|. Insets show the group velocity Vg — duj/dq as a 
function of lj. 



From this fit we can find both the phonon frequency 1) we can write instead of Eq. (10| 



Wq and the phonon line width Aw. The results for ui^ are 
shown on Fig. [7]for three values of /x and q || (100). For 
sufficiently small values of wave vector q we see a nice 
linear dispersion curve ujq = vq, with the sound velocity 
V given by Eq. ([9). It is independent of the q direction 
(i.e. the sound velocity is isotropic). With increase of q, 
the frequency uiq shows a pronounced negative dispersion 
and approaches the maximum frequency Wmax where the 
dependence ujg saturates. In this q region we observed 
a weak anisotropy of the dispersion curves for /x = 1. 
At smaller values of /i the dependence Wq is isotropic. 
Since Wmax oc ^JJi, the vertical axis on Fig. [7] scales ap- 
proximately as ^JJl and the horizontal axis scales as /i^/"* 

(sound velocity v oc ^/E oc /i^^'*, and Qmax ~ i^max/f ). 

The strong negative dispersion in ujq for big q values 
can be explained by term crossing effect due to the cou- 
pling of phonons to quasilocal vibrations near frequency 
'^maxi corresponding to the sharp maximum in the DOS 
g{ui) (see Fig.|4]). The dip in the participation ratio P{uj) 
for ji = 0.001, /X = 0.01 and /i = 0.1 at w w„iax (see 
Fig. [5]) evidences in favor of this idea. The vibrations 
inside the dip correspond to frequencies near Wmax and 
have smaller participation ratio than the others. There- 
fore they can be referred to as quasilocal vibrations. In 
the following we will see that this strong scattering is 
also responsible for the deep minimum in the diffusivity 
D[lS) at w sa Wniax (see Fig. 17 1. Finally let us note that 



observed negative dispersion has nothing to do with the 
negative dispersion of transverse phonons in crystals near 
the Brillouin zone boundary. 

The negative dispersion in Uq is responsible also for 
the pronounced rise of the phonon DOS above the 
dependence, given by Eq. ( |lOl ). It is clearly seen on the 
Fig. |4] Indeed, taking the dispersion into account and 
disregarding weak anisotropy (taking place only for /i = 



5ph(w) 



1 q\^) 

27r2 Vg{u!) 



(15) 



Here Vg — du j dq is the group velocity shown in Insets on 
Fig. [7] Taking for qiuS) and Wg(a;) the data from Fig. [t] we 
obtain the points (filled and open diamonds) shown on 
Fig. |4] Since they perfectly coincide with numerical data 
for g{y)) below Wmaxj we conclude that aHthe excitations 
in the phonon gap belong to phonons (with nonlinear 
dispersion at higher values of q) . 
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FIG. 8: The phonon line width Aa; as a function of a; for 
different /i in cubic sample with A'^ = 50^ (one realization). 
Different symbols correspond to different q directions. O for 
q II (100), A for q || (110), □ for q || (111). Filled and open 
symbols refer to excitations below and above the loffe-Regel 
crossover frequency ojir correspondingly (see text for details) . 

The phonon line width Aw can be also found from fits 
similar to those shown on Fig. [6] It is related to the 
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phonon life time r = l/2Aa;. The factor 2 takes into ac- 
count that Aw corresponds to decay of the amphtude of 
the vibration. The results are shown on Fig. |8] As follows 
from this figure, Auj oc cj'' and does not depend on the 
direction of q. We think that this frequency dependence 
is not due to Rayleigh scattering of phonons on a static 
disorder. In such a case Aw would be proportional to q^. 
Due to nonlinear dispersion in Wq, these dependencies do 
not correspond to each other. More likely, the phonon 
line width is due to strong resonant scattering of phonons 
by quasilocal vibrations responsible for the sharp peak in 
the DOS, similar to those introduced irP. The deep min- 
imum in the diffusivity D{uj) around frequency Wmax also 
supports this idea (see Fig. [T7|) . We hope to investigate 
this important question in future work. 

With known value of Aw, the phonon mean free path 
Z(w) can be calculated as follows 



/(w) 



^ 2Ao 



(16) 



The phonons are well defined excitations if their mean 
free path /(w) exceeds the phonon wave length A = 27r/g 
(lofFe-Regel criterium for phonons). As we will see in 
the next Section, phonons transform to diffusons when 
w A/2. We will call the corresponding crossover fre- 
quency as wiR. Fig. |9] shows the ratio l{u>)/X as a func- 
tion of w for several values of fi and different directions 
of the wave vector q. The boundary between filled and 
open symbols (the full horizontal line) corresponds to fre- 
quency wiR. Thus filled and open symbols on Figs.|4][7j 
|8j|9] belong to phonons with frequencies below and above 
the loffe-Regel crossover frequency correspondingly. 
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FIG. 9: The ratio l{uj)/X as a function of uj for different 
fi. Different symbols correspond to different q directions 
as explained on Fig. [s] The full horizontal line (separating 
filled and open symbols) corresponds to loffe-Regel crossover 
l{uj) = A/2. 

Usually in glasses the loffe-Regel crossover frequency 
wiR is correlated with position of the boson peak Wf,, 
gg^HHm Qjr^^ references therein. It is the frequency where 
the reduced DOS g{u})/uj^ has a maximum. We also have 
a rather sharp boson peak in our disordered lattice^^. 



As is evident from Fig. [4j it comes from vibrations in the 
low frequency range of the dynamical matrix M = AA^ . 
These vibrations are pushed out to higher frequencies 
when we add to the dynamical matrix AA^ a crystalline 
part (iMq. Just these vibrations, forming an additional 
density of states near frequency Wmax, are responsible for 
the boson peak in our case. As follows from Fig. |4] the 
left side of the boson peak is constructed from phonons 
having negative dispersion Wq. Similar results were ob- 
tained recently for Lennard- Jones glassej^* **^ * "^^ *. It is not 
inconceivable that the boson peak in some other glasses 
can have a similar structure. 

The frequencies Wmax, wir, and Wb are collected in Ta- 
ble H] for different /i. As we can see from the table, wjr 
is close to the frequency w^ax and to the position of the 
boson peak uJb- Above wir phonons cease to exist as 
well defined excitations. They are smoothly transformed 
to diffusons which we will consider in the next Section. 
The relative number of phonons in the lattice can be es- 
timated as follows 



g{uj)duj. 



(17) 



These values are also given in the Table |I] We see that 
for all investigated values of fi the relative number of 
phonons in the lattice is small. It is in agreement with 
similar estimates for amorphous silicorP^. 









torn. 


iVph 


1 


2.5 


2.4 


2.2* 


0.12 


0.1 


0.78 


0.74 


0.62 


0.027 


0.01 


0.23 


0.23 


0.19 


0.0066 


0.001 


0.072 


0.07 







TABLE I: The frequency of maximum in DOS cJmax, the fre- 
quency of the loffe-Regel crossover ljir and the boson peak 
frequency uji, for various jj,. Star * means that cjir was found 
for q II (100) direction. A^ph is a relative number of phonons 
in the lattice. 



IV. DIFFUSONS 

In this section we are going to consider properties of 
diffusons. As is well known, the diffusion phenomenon 
usually takes place for physical quantities which are con- 
served. In a free closed mechanical system we have two 
integrals of motion, momentum and energy. 



A. Diffusion of momentum 

First let us consider diffusion of momentum. Usually 
the diffusion of momentum is related to viscosity in the 
system. All particle masses being equal (m^ = 1), the 
diffusion of momentum is equivalent to the diffusion of 
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FIG. 10: The displacement structure factor 5(q,a;), Eq. (symbols) for ^ — Q and for three different frequencies. The sample 
size IS iV = 50^ The averaging is performed over 300 realizations. Different symbols correspond to different q directions. O 
for q 11 (100), A for q || (110), □ for q |j (111). Full lines correspond to the structure factor 5'rw(q,'^) of the random walk on 
the lattice given by Eq. (19 1 with _Drw = 0.7. Dashed line corresponds to the limit g ^ 1 (see Eq. (22|). 



particle displacements. It is because in our system the 
position of the center of inertia is conserved and we can 
put it at the origin of the coordinate system. Then the 
sum of all particle displacements vanishes 



^M,(i) =0, 



(18) 



i.e. it is an integral of motion. The diffusion of displace- 
ments in this case looks like a diffusion of "particles" in a 
lattice vifhere the total number of particles is conserved. 

By analogy with diffusion of "particles" the informa- 
tion about diffusivity of displacements is absorbed in the 
displacement structure factor S{(\,uj) ( |Tl| ). We remind 
that to calculate this structure factor we ascribed at the 
initial moment t — the random displacements to all the 
particles with Gaussian distribution (with zero mean and 
unit variance) and velocities equal to zero. So the condi- 
tion ( 18 ) at i = was satisfied. Therefore let us analyze 
now this structure factor in the diffuson frequency range. 

Consider first the case of /i = when phonons are ab- 



sent and only diffusons are present in the lattice. Fig. 10 
shows the structure factor S'(q, w) as a function of wave 
vector q for three different directions in q space (symbols) 
and for three different frequencies ui. Let us compare this 
displacement structure factor with structure factor of the 
random walk S'rw(q, w) on the lattice. 

As was shown irP^l for the case of the random walk on 
a lattice, S'rw(q, w) is given by expression 

2r(q) 



r2(q) 



(19) 



It is a Lorentzian, with a width r(q) given by 
r(q)- AwQ'(q), 



(20) 

where i'rw is a diffusion constant of the random walk. In 
a simple cubic lattice (with lattice constant — \) the 
fimction Q{(\) reads 



(3(q) = 2Wsin 



sm' 



2 



sm 



(21) 



For small values of q ^ 1, Q(q) = q and in the contin- 
uum limit we have the well known result for the diffusion 
structure factor 



S'rw(q,t^) = 



(22) 



Let us note that the structure factor ( 19 1 has a maximum 
at q values obeying the condition 



r(q) - Z?rwg'(q). 



(23) 



We can specify it as a dispersion law for diffusons. The 
width of the maximum is r(q). For g ^ 1, r(q) = D^y^q^ . 




FIG. 11; The correlation function C(r, cj) for /i = and six 
frequencies u (0.14, 0.31, 0.49, 0.66, 0.84, 1.01) for sample 
with A'' = 50^ particles averaged over 300 realizations. The 



full lines are our numerical results obtained from Eq. (111. 
Each line starts from r = rmin which is about 2.5 interatomic 
distances (marked by arrows). The dashed line corresponds 
to Eq. p6| with = 0.7. 



A comparison of the displacement structure factor 



5(q, oj), (111, and the structure factor of the random 
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FIG. 12: The normalized structure factor ^^(q, uj) as a function of q for some direction in q space and for each frequency uj for 
various values of fi (0, 0.01, 0.1, 1). The sample size is N ^50^. The averaging is performed over 100 realizations. Left sides 
of all plots are for q || (111), right sides are for q || (100). White horizontal dashes show the lofle-Regel crossover frequency 
ojiR. For fJ, = 1 the frequency uuib. is slightly different for different q directions. Black full line corresponds to Eq. (231 for the 
random walk on a simple cubic lattice with diffusion constant Drw = 0.7. 



walk S'rw(q, '^), (19), is shown on Fig. 10 One fitting 
parameter was the diffusion coefficient in Eq. (20 1. 



From comparison of these data we obtained sa 0.7. 
It means that the diffusion coefficient of particle displace- 
ments Du « 0.7 (see Section|v]). Another fitting parame- 
ter was a height h{uj) of the random walk structure factor 
in the maximum. According to Eq. ( 19 1, in the maximum 
r(q) — CO and h{uj) = l/oj, but to fit the data points on 



Fig. 10 we used slightly higher values of h{uj) 



The small difference between /i(a;) and 1/w can be ex- 
plained by different frequency dependencies of the den- 
sity of states g{uj) for vibrations and for the random walk 
(following from the sum rule similar to Eq. ([l3|). As 
we can see from the figure, for the investigated frequen- 
cies the fit is perfect. With increasing frequency above 
w « 2 — 3, the fitting becomes more and more poor since 
we approach the localization threshold at wioc ~ 5.5 ±0.5 
(see below) which is not described well by a simple model 
of Markovian random walk. 

Now let us consider a behavior of a correlation func- 
tion. The correlation function of particle displacements 
at some frequency uj, expressed through eigenvectors 
er{u)) of the dynamical matrix M, reads 



C(r,w) = ^ er'+r(a;)er'(a;). 



(24) 



It is a Fourier transform of the displacement structure 
factor (111 



^ J 5(q,a.)e^'i"-dq. 



(25) 



Let us compare this correlation function with correla- 
tion function of the random walk. For distances bigger 
than the period of the lattice (oq = 1) we can make use 
of the limit of small q 1 and integrate Eq. ( 22 ) for the 



random walk structure factor taken in approximation of 
continuous media. As a result, we derive 



exp 



Crw(r, w) = 



2t:'^ rD„ 



2Dr 



(26) 



Fig .|ll| shows a good agreement of our correlation func- 
tion ( 25 1 with the correlation function of the random 
walk ( 26 1 . For all investigated frequencies the numeri- 



cal data collapse together and become indistinguishable 
from the theoretical prediction ( 26 1 . We can see also on 



this figure the anticorrelation phenomenon (the region of 
negative values of the correlation function). As follows 
from Eq. ( 26 ) , the correlation function of the random 



walk changes its sign for the first time at 



UJ 



TT 

2' 



(27) 



It is also in a good agreement with our numerical results. 
Therefore we can call a corresponding value of r found 
from Eq. (27) as a radius of diffuson. It is a typical size 



of the regions vibrating with frequency uj and having the 
same sign of all particle displacements. 

Now let us analyze the displacement structure factor 
5(q, uj) for 0. For better visual effect we will show a 
map of the function S'(q, w) on the plane (w, q) for differ- 
ent directions in q space. To do that, for each frequency 
UJ we have found the maximum S'(q, uj) as a function of 
q along some directions in q space. Then we normalized 
function S'(q, uj) along this line cj=const to the magni- 
tude of this maximum. 



The results are shown on Fig. 12 for four different val- 
ues of fj, and two directions in q space. The white color 
corresponds to the maximum when normalized structure 
factor S'„(q, w) — 1 while the black color to the case 
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where S'„(q, w) = 0. For /x 7^ we can see clearly two 
types of excitations in the lattice. At low enough fre- 
quencies, below wiR, we see phonons with well defined 
dispersion law cjq, the same as in the previous Section. 
At the loffc-Regel crossover frequency wir , the structure 
factor strongly broadens and phonon dispersion line dis- 
appears. Above the displacement structure factor co- 
incides well with the structure factor for /i = case shown 



200 



on Fig. 12 1, which corresponds to diffusons. The maxi- 
mum of the normalized structure factor S'„(q, cj) (white 
regions) agrees well with Eq. ( 23 ) giving the maximum of 
the random walk structure factor S'rw(q, (black line). 
Deviations from 5'rw(q, w) take place at high frequencies 
near the localization threshold. 
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FIG. 13: The same normalized structme factor ^^(q, a;) as 
on Fig. 12 but in q space in plane q^qy {qz ~ 0) for cj = 0.5. 
The left picture corresponds to fi = (a) and the right to 
fi = 0.1 (b). 



Finally, to compare phonon and diffuson structure fac- 
tors, a cross section of the structure factor S'„(q, w) in 
q space for = and frequency lu = 0.5 is shown on 
Fig. 13 for /i = and /x = 0.1. At the left side (a) of 
this figure we see the structure factor of diffuson. On 
the right side we see the structure factor of phonon (b). 
As compared with phonon structure factor, the diffuson 
structure factor is much more broadened. 



B. Diffusion of energy 

Now let us consider the diffusion of energy. The diffu- 
sion of energy is different from diffusion of particle dis- 
placements (see Section |v]). The first approach to cal- 
culate the diffusivity of energy D(uj) for vibrations with 
frequency w is a direct numerical solution of Newton's 
equations. For that we have used the Runge-Kutta-4 
method with time step At = 0.01 applied to a cubic 
sample with N — L x L x L particles (lattice constant 
ttQ — 1) and with free boundary conditions along the x di- 
rection. Along other two directions we take the periodic 
boundary conditions. 

Assuming zero initial conditions for displacements and 
velocities of all the particles, let us apply external forces 
with frequency uj and random phases ifi to the particles 




1000 



FIG. 14: The dependence of R^{t) in the case of /i = for 
one sample with A'' = 100 x 100 x 100 particles and 14 dif- 
ferent frequencies uj — 0.5, 1, 1.5, . . . , 7 (from top to bottom). 
The numbers indicate integer frequencies. The slope of each 
line corresponds to each black dot in Fig. [15] Two points at 
(jj — 2 and uj = 6 correspond to two distributions of energy 
E{x, t) over the sample for delocalized and localized modes 
correspondingly. They are shown on Fig. 16 (see below). 



in the central layer a; = of our sample 



f-''\t) =sin(a;t-|-(^,)exp 



2r2 J 



(28) 



where utT ^ 1. The right and the left sides of the sample 
have coordinates a;r,i — ±L/2. In such a way we excite vi- 
brations with frequencies near frequency uj distributed in 
a small frequency interval (w — 1 /T, uj + l/T). In calcula- 
tions we used T = 5 for all frequencies uj. We started our 
calculations at time to = —5T when the external force is 
still negligible. 

After applying the force to the central layer a; = 0, 
vibrations will spread to the left and to the right ends of 
the sample. The average squared distance to the energy 
diffusion front we define as usual 



1 

E7r 



N 
i=l 



1 



L/2 



E,. 



x^E{x,t)dx. (29) 



-L/2 



Here Xi is the x coordinate of the i-th particle, Ei{t) is the 
energy of i-th particle and sum is taken over all particles 
in the sample. E'tot = J2i Ei{t) is the total energy of the 
system. It is independent of time after the external force 
/™*(t) becomes negligibly small (i.e. for t > 5T). 

The energy of «-th particle Ei {t) we define as a sum of 
the kinetic energy and a half of the potential energy of 
connected bonds (m; — 1) 



2 



(30) 



Here Vi{t) = Ui{t) is a particle velocity. Summation over 
all particles in Eq. ( 29 ) we can divide in two steps. First 
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we sum over all particles in the layer x and then we sum 
over all layers. Let E{x,t) be a total energy confined to 
the layer x at time t. Having in mind that in our case we 
have lattice constant oq = 1 and sample size i 3> 1, we 
can change summation over different layers to integration 
over coordinate x for times where R{t) ^ 1. 



— 10^ Thouless 

— 14^ Thouless 
20^ Thouless 

• 100^ Newton 




:3 







->-to = 


2 






L — 




- 




-*- ai = 


6 - 






V ^ ^ 


900 




/A 







-50 



50 



FIG. 16: Black points (diamonds and triangles) show the dis- 
tribution of energy E{x, t) contained in the layer x as a func- 
tion of X for two different frequencies uj = 2 and tj = 6 at times 
t = 234 and t = 900, respectively, calculated numerically with 
Newton method. Full lines are theoretical predictions for de- 
localized (diffusive) and locahzed modes given by Eqs. (32 
33 1 with ^ 166 and « 22 correspondingly. 



FIG. 15: The dependence of diffusivity D{uj) on uj for /i = 0. 
Black dots are calcul ated by the direc t so lution of Newton's 
equations from Eqs. (29|[3l| and Fig. fT4|for A*' = 100^ par- 
ticles (one realization) 



FuU lines for iV = 10^14^20'' are 



calculated using formula of Edwards and Thouless ( 36 1 with 
c = 1 (see below). Averaging for lines is performed over fre- 
quencies in the small interval {uj — 5uj, uj + 5uj) with 5uj = 0.25 
and over several thousands realizations. 

We will apply this method to the case of /i = (i.e. 
for the lattice without phonons). The results are shown 
on Fig. 14 As we can see from the figure for small and 



middle frequencies, R {t) oc t. Therefore for these fre- 
quencies vibrations indeed spread along the x axis by 
means of diffusion. The slope of the lines decreases with 
frequency w. For calculating the slope, we take the time 
interval At where, on the one hand t > 5r, and on the 
other hand, R < L/2. 

From the slope of R^ (t) we can calculate the diffusivity 
of modes D{uj) using one dimensional formula 



R^{t) = 2D{uj)t. 



(31) 



This diffusivity is shown by black dots on Fig. [T5| At 
small frequencies it is approximately constant, then it 
decreases with frequency approaching zero at the local- 
ization threshold, wioc « 5.5 ± 0.5. At higher frequen- 
cies above wioc the dependence R^{t) saturates with in- 
creasing t. This indicates localization of the vibrational 
modes. 

The difference between delocalized and localized modes 
is clearly seen if we examine the dependence E{x, t) as a 
function of coordinate x at some moment t for two differ- 
ent frequencies below and above the localization thresh- 
old. These two points for investigation are shown on 
Fig. [l4j Black diamond corresponds to delocalized mode 
with frequency uj = 2 and has coordinates t = 234 and 
R^ = 166. The distribution of energy E{x,t) over the 



sample calculated numerically at this moment is shown 
by black diamonds on Fig. [THj The data are perfectly 
fitted by solid line drawn according to the solution of 
diffusion equation in Id case 



E{x,t) = 



Et. 



V2ttR^ 



exp 



x^ \ 
'2i?2 J 



(32) 



with value of R^ ~ 166. 

Black triangle on Fig. [T4|corresponds to localized mode 
with frequency w = 6 and has coordinates t — 900 and 
R^ = 22. The distribution of energy E{x, t) over the sam- 
ple calculated numerically at this moment is shown by 
black triangles on Fig. [16] This distribution is drastically 
different from the previous case. For localized modes we 
expect the usual exponential decay 



E{x,t) 



Et, 



/2R 



exp 



V2\a 
R 



(33) 



The fit of the numerical data with this function and 
R^ = 22 is shown on Fig. 16 The fit is perfect except 
for the central point at x — which lies noticeably above 
prediction of Eq. (33 1. The coefficients in Eqs. (32 



33) 



were taken to satisfy the obvious rules 

oo oo 

/ E{x,t)dx = Etot, -^r- I x'^E{x,t)dx = R^. 
J Etot J 

— oo —CO 

(34) 

To find the diffusivity D{uj) for /i ^ 0, the method 
of numerical solution of Newton's equations is not ap- 
propriate, because in this case we have phonons in the 
lattice with long mean free paths. Correspondingly sam- 
ples with much bigger sizes are necessary to use this ap- 
proach. Therefore for /i 7^ we used a second approach. 
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In this approach, the diffusivity D{uji) at eigenfrequency 
LOi was calculated by means of the formula of Edwards 
and ThoulesPI 



For fj, = the results for D{uj) are shown on Fig. 15 for 



(35) 



where L is the length of the sample and Awi is sensi- 
tivity of the eigenfrequency iOi to a twist of boundary 
conditions. More precisely, we defined the diffusivity as 
follows: 



DH = clim — (lA^HI) 



(36) 



where (p is the angle of twisting, and c is some constant of 
the order of unity. It will be determined from comparison 
with the Newton method. The averaging in Eq. |36] is 
performed over frequencies uj in the small interval (w — 
duj,uj + 6uj) with Slu = 0.25 and/or over several thousands 
realizations. 




a> 



FIG. 17: The diffusivity D{lu) for various fi (0, 0.01, 0.1, 1) for 
sample with A'^ = 14"^ (crosses). The diffusivity was calculated 



using formula of Edwards and Thouless ( |36[ ) with c — I and 
averaged over two thousand realizations. The arrows indicate 
frequencies cJmax in the DOS g{ijj) for correspondin g va lues of 
jj. Open symbols correspond to phonon diffusivity ( 39 1 below 
the loffe-Regel crossover frequency ujir. 

The symmetric real matrix M was defined as usual ^ 
with periodic boundary conditions. The twisting of the 
matrix M by angle ip gives a new Hermitian matrix M' 
obtained as follows. For bonds between the left (l) and 
the right (r) boundaries of our cubic sample 

A4 = Mir exp(i^), = Mri cxp{-iip). (37) 

For all other bonds Mj^, = Aljk- So Auji is the difference 
between i-th eigenfrequencies of matrices M and M' 



Aw,; 



- W,- 



(38) 



Twisting of boundary conditions was performed for x 
direction only. For others two directions the periodic 
boundary conditions were used. 



three different cubic samples (full lines). We compared 
these results with numerical solution of Newton equations 
for = (black dots) and get for the constant c « 1. 
Then we used this c value for /i ^ 0. The results are 
shown on Fig. [TT) For /i 7^ we see clearly two different 
frequency regions in the function Diuj). 

At low frequencies, diffusivity increases with decreas- 
ing of UJ. This range corresponds to the phonons. Indeed, 
the diffusivity of phonons D{uj) can be calculated as fol- 
lows 



1 



D{uj) — -l{uj)vg{uj). 



(39) 



Open symbols on Fig. 17 show contribution calculated 
from this equation (just below loffe-Regel threshold). We 
see a good agreement with Edwards and Thouless for- 
mula. After a deep minimum at frequency uj » Wmax 
the diffusivity D{uj) saturates at a constant level coin- 
ciding with D{uj) for ^ — Q. The diffusivity in this range 
corresponds to diffusons . Sim ilar behavior of D{ijj) was 
found in jammed systemj^lllll. The deep minimum in the 
diffusivity at a; w t^Jmax corresponds to strong scattering 
of phonons by the quasilocal vibrations near the sharp 
peaks in the DOS g{uj) (see Fig. [I]). 



V. DISCUSSION 

We have developed a random matrix approach to de- 
scribe vibrations in strongly disordered systems, which 
have properties similar to what one observes in granular 
matter at the jamming transition point, in jammed sys- 
tems and, finally, in real glasses. This approach has one 
important advantage in comparison to other models. It 
describes mechanical systems which are always stable in- 
dependently o f the de gree of disorder. Previous random 
matrix models'^^'^SlIl] suffer from an inherent mechanical 
instability that occurs at some critical amount of disor- 
der. As a result, they are limited by consideration of 
"relatively weak" or "moderate" disorder. 

We take the dynamical matrix in the form M — AA^ + 
fj-Mg. Here ^ is a random matrix N x N built on a simple 
cubic lattice with N particles and interaction between 
nearest neighbors only. The only non zero non-diagonal 
matrix elements Aij between the nearest neighbors are 
taken as independent random numbers from Gaussian 
distribution with zero mean (Aij) — and unit variance 
{A'^j') = 1. The variance controls the degree of disorder 
in the lattice. To ensure the translational invariance the 
diagonal elements are calculated as a minus sum of non- 
diagonal elements An = —J^j^^i ^ji- ^0 is a simple 
crystalline dynamical matrix with unit springs between 
the nearest neighbors. 

If the first term AA^ is responsible for the disorder in 
the system, the second term fiM^ describes the ordered 
part of the Hamiltonian. The parameter fi controls the 
relative amplitude of this part and the rigidity of the 
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lattice. It can vary in the interval < < c«, changing 
the rigidity and relative amount of disorder. In this paper 
we have mainly considered the case of strong disorder 
when < /i < 1 and fluctuating part of the dynamical 
matrix is bigger then the ordered part. In this case the 
Young modulus of the lattice E cx ^/JI. Such form of the 
dynamical matrix guarantees the mechanical stability of 
the system for any positive value of fj,. 

We have found that the delocalized vibrational exci- 
tations in this disordered lattice are of two types. At 
low frequencies below the loffc-Regel crossover, w < wir, 
they are the usual phonons (plane waves) which can be 
characterized by frequency uj and wave vector q. How- 
ever, with increasing of ui, due to the disorder-induced 
scattering, the phonon line width Aw increases rapidly 
as Auj (X Lu* and at some frequency uj « wir the phonon 
mean free path I becomes of the order of the wave length 
A. Though this crossover is not sharp and has no critical 
behavior at a; = i^ir, the structure of the eigenmodes at 
higher frequencies quite soon become very different from 
the plane waves. 

As a result, at higher frequencies the original notion of 
phonons is lost and delocalized vibrational modes have a 
diffusive nature. They are similar to dijfusons introduced 
by Allen and Feldman, et alW^. The diffusons again can 
be characterized by frequency w, but have no well defined 
wave vector q. Above lo « wir the structure factor of 
particle displacements S{q, uj) becomes very similar to 
the structure factor S'i.„(q, w) of a random walk on the 
lattice. The former has a broad maximum as a function 
of q at q ~ ^^JJ|Du^ where is a diffusion coefficient 
of the particle displacements. 

The displacement structure factor S((i, t) in the dif- 
fuson range, for small q ^ 1/ao, decays as following, 
S{q,t) (X exp{—Duq^t). As a result the vibrational line 
width T{q) = D^q^ . Such quadratic dependence of T(q) 
was found in many glasses in the experiments on inelas- 
tic x-ray scattering, see for examplefS^ES and references 
therein. It was also found in molecular dynamic simula- 
tion of amorphous silicorf^. However in these and other 
papers this line width was attributed to phonons with- 
out discussion of its physical origin. We guess that the 
observed q^ dependence of T{q) has nothing to do with 
phonons and is in fact related to diffusons. However, 
a more detailed investigation is necessary for a definite 
conclusion. 

The crossover between phonons and diffusons takes 
place at the loffe-Regel crossover frequency wir which 
is close to the position of the boson peak. Since for 
phonons Aw oc and for diffusons T{q) — Duq^, there 
should exist a crossover from uj* to dependence of the 
line width. Such a crossover was indeed found recently 
in inelastic x-ray scattering in lithium diborate glass^, 
densified vitreous silica^^, vitreous siliccP^^, glassy sor- 
bitol^ and glycerol glas^^. The crossover frequency was 
found to be close to the BP position. 

As a result, if our guess is true, we can calculate 
the diffusion coefficient of particle displacements, = 



r{q)/q'^, from the experimental line width r{q) in the 
range, where it is proportional to q^. Taking into ac- 
count that ~ ^o/''' where ap is the lattice constant 
and r is an average time for a jump, we come to the or- 
der of the value estimate Du ~ Imm^/sec for ao ~ 2A 
and T « 0.4 x 10^ ^'^ sec. Let us compare this value with 
experimental data. 

In the papei'^^ it was found that in vitreous silica 
hT/{hujf = O.OTmeV^i for q >2 nm^i. Taking the 
sound velocity = 5250 msec^^ for q = 2 nm~^ we 
get for diffusion coefficient D„ = 1.3mm^/sec. Let us 
compare this value with the diffusivity of energy D{uj) 
for small uj in the same glass. We expect that both coeffi- 
cients should be of the same order of magnitude. The dif- 
fusivity of energy D{uj) in vitreous silica was calculated in 
the papejl^Sl. It was obtained that D{0) = 1.4 mm^ /sec. 
A close estimate D{0) ~ l.lmm^/sec was given irPSl. 
As one can see the agreement between Du and D{0) is 
unexpectedly good. In glycerol glas^^ we found the dif- 
fusivity about factor of two smaller, D„ = 0.46mm^/sec. 
For amorphous silicon from molecular dynamic calcula- 
tioni^we get D„; = 3.2 mm^ /sec for longitudinal vibra- 
tions, and Dut = 1.2mm^/sec for transverse vibrations. 
For the diffusivity of energy we have in this glass the 
estimate^^ ^(0) — 0.6nim^/sec. 

Since wir c>c ^/Jl, we can vary the loffe-Regel crossover 
frequency and, therefore, the relative number of phonons 
TVph in the system, changing the parameter /i. It is zero 
when fi = and there are no phonons in the lattice. 
In this case all delocalized vibrations are diffusons. If 
< /i ^ 1 we have phonons, but their relative number 
is small. One can show that in this case Npi^ oc fi^^'^. In 
the opposite case, /i ^ 1, the disorder is relatively small 
and nearly all vibrations in the lattice are well defined 
plane waves, i.e. phonons. 

In amorphous silicon the relative number of phonons 
(plane waves) was estimated to be only 4% from all of 
the vibrational modes in the systerrpSl. The estimates 
show that in our model we have such a small amount of 
propagating modes, as in a-Si, for 0.1. In the silica 
glass we can estimate the relative number of phonons 
from the dataP^'. Taking into account that loffe-Regel 
crossover frequency in amorphous silica was estimated to 
beP fiR = ITHz, and integrating density of states" up 
to this frequency we come to the relative number A^ph — 
0.002 ± 0.0005. As a result in the typical glass such as 
amorphous silica only 0.2 % of all modes are phonons. As 
follows from Table [l] it corresponds to very small values 
of /i < 0.01. It means that small amount of phonons in 
disordered systems is a signature of strong disorder. 

Usually the phenomenon of diffusion takes place for 
conserved quantities. In our system we have two inte- 
grals of motion. They are the momentum and the energy 
of the lattice. Therefore, first of all, one has to discrimi- 
nate the diffusion of particle momentums (or particle dis- 
placements) from the diffusion of energy. Conservation 
of displacement is related to conservation of the center 
of inertia in the system. As a result, the diffusion of par- 
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tide displacements has the same diffusion coefficient as 
the diffusion of particle momentums. 

The diffusion coefficient of displacements/momentums 
D^iy is hidden in the displacement structure factor 



5(q, w) (11). Comparing this structure factor with the 
structure factor of the random walk on the lattice, we 
found that for the case of /i = the diffusion coefficient 



D 



u/v 



0.7. We can check that it is indeed the dif- 



fusion coefficient of particle displacements/momentums 
in a similar way we used for finding the diffusivity of 



energy D{uj) in Section IV B 



Let us consider a cubic random lattice L x L x L with 
H = and unit masses nii = 1 with periodic boundary 
conditions. At initial moment t = let us displace all 
particles in a thin layer around the central layer (with 
coordinate a; = 0) according to Gaussian distribution 



(40) 



Here the thickness of the layer xq should be small enough 
in comparison to the sample size L, i.e. xq ^ L/2. Initial 
velocities u{0) of all the particles are equal to zero. 

After initial displacements in the thin central layer, 
the particle displacements will diffuse to the left and to 
the right ends of the sample. Solving numerically the 
Newton equations, we find the average squared distance 
to the displacement diffusion front, similar to Eq. (29) 



Rlit) 



—'^xfuiit), Mtot = y^u»(t). 



(41) 



Since the center of inertia does not move, the total dis- 
placement of all particles Utat is independent of time and 
equal to the total displacement a,t t — 0. 

From the slope of R^{t) we can calculate the diffusion 
coefficient of the displacements Z)„ as follows 



Rlit) = 2Dut 



(42) 



similar to Eq. (31 ) 



In the same way we can calculate the diffusion of mo- 
mentum. For that at the moment t = initial displace- 
ments of all the particles we put equal to zero. However 
initial velocities v — u(0) in the thin central layer we take 
distributed similar to Eq. ( 40 ) 



v{x) — v^e 



(43) 



Then, as in the previous case, solving numerically the 
Newton equations we find 



Rlit) 



1 

Vtot 



xlvi{t), Vtot 



E 



v,it). (44) 



Since the total momentum is conserved, Wtot is also in- 
dependent of time and equal to its initial value at i = 0. 
From the slope of Rl{t) we can calculate the diffusion 
coefficient of the momentum using one dimensional 
equation 



Rlit) = 2D J 



(45) 



similar to Eq. ( 42 ) 



In both cases we have obtained for diffusion coefficients 
Du and £)„ the same value as was derived from the struc- 
ture factor, Du ~ « = 0.7. It confirms our state- 
ment that the displacement structure factor S'(q, w) gives 
us the information about diffusion of particle displace- 
ments (or momentums). The diffusion of momentum is 
usually related to viscosity rj of the medium. Therefore 
in the case of /i = our lattice has no rigidity but has a 
finite value of viscosity. 

In disordered lattices the diffusion of energy is differ- 
ent from the diffusion of particle displacements (momen- 
tums). In the harmonic approximation the eigenmodes 
with different frequencies do not interact with each other. 
Therefore the energy cannot be transferred from one 
eigenmode to other eigenmodes. It means that energy 
of every eigenmode Eiuji) is conserved (with time). The 
total energy i?tot is just a sum of these eigenmode con- 
tributions 



(46) 



As a result, instead of one integral of motion (the total 
energy i?tot), in a scalar harmonic system with N parti- 
cles we have N integrals of motion E{uji). And for each 
frequency we have its own unique energy diffusivity 
D{uji). At this point our model decidedly confirms the 
physical picture suggested in papers^^ ■^^ for amorphous 
silicon. We believe that it can be applied to some other 
glasses as well. 

Usually this diffusivity is hidden in a displace- 
ment/momentum structure factor of the 4-th order. 
However, we calculated the diffusivity of energy -D(w) in 
a different way using two different approaches as it was 



discussed in Section IV B The first approach is based 



on the direct solution of Newton equations. In the sec- 
ond approach we calculated the diffusivity using Edwards 
and Thouless formula^. Both approaches give the same 
result. 

In the first approach we used a short external force 
pulse At exciting vibrations in a small space region of the 
lattice and in a small frequency interval Aw 1/ At near 
frequency uj. Then on a time scale t 3> At the energy 
diffused through the lattice. Using Newton equations 
of motion we calculated this diffusion directly. It was 
supposed that the interval Auj is much bigger than the 
interlevel spacing Suj and therefore the former consists 
of many eigenmodes. In the thermodynamic limit 6uj oc 
1/N — )■ if N oo. Therefore in an infinite system 
we can take the interval Auj arbitrary small. The energy 
diffusion coefficient D{ll!) in this case is a function of 
frequency uj. Approaching the localization threshold ujioc 
the diffusivity I?(w) should go to zero. 

We applied this method for ^ = 0, when there are no 
phonons in the lattice. In this case we obtained for diffu- 
sivity at zero frequency 13(0) ~ 0.4, i.e. the value about 
factor of two smaller than for diffusivity of displacements, 
Du- However this approach is rather difficult to imple- 
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ment for computer simulations in the case when fi 0. 
In this case we have phonons in the lattice with long 
mean free paths. And samples with much bigger sizes 
are necessary. 

Therefore, to calculate the diffusivity D{uj) for arbi- 
trary value of /i (including the case of /i = 0), we used 
another approach. In this approach Edwards and Thou- 
less formula^'^, D{uJi) = cL^jAwij, was used. It relates 
the diffusivity D{uji) with shift of the eigenfrequencies 
Auji due to change of the boundary conditions in one di- 
rection. The proportionality coefRcient c we found from 
the comparison with the Newton method for /i = 0. In 
this case both methods result in the same frequency de- 
pendence of D{uj). 

The diffusivity of vibrational modes D{ll)) in disordered 
lattices is a very important quantity. It determines the 
thermal conduct ivity^SI 

CO 

k{T) j dujg{ijj)D{uj)C{uj,T). (47) 



Here 5(1^) is density of states and C(aj, T) is specific heat 
of harmonic oscillator 



(48) 



Localized modes have D{uji) = and make no contribu- 
tion to x{T). 

If functions g{uj) and D(uj) are approximately constant 
in some frequency interval (the case that we have, for ex- 
ample, in our picture for oj > ujm), then we find from 
Eq. 47 that approximately >f(r) cx T in the correspond- 



ing temperature range. It explains a quasi-linear temper- 
ature dependence of the thermal conductivity above the 
plateau observed in glassed. With increasing frequency 
the functions g{uj) and Diuj) finally drop to zero and 
thermal conductivity saturates at some constant level in- 
dependent of temperature. Thus the conception of diffu- 
sons gives clear explanation for the temperature depen- 
dence of the thermal conductivity of glasses and other 
disordered systems. 

Summarizing, using a stable random matrix approach 
we have presented a consequent theory of vibrational 
properties in strongly disordered systems. In these sys- 
tems a relative amount of phonons is small and almost all 
delocalized vibrations are diffusons. The diffusons play 
an important role and are responsible for the transport 
properties of glasses at higher temperatures. Presum- 
ably they are also accounted for the mysterious de- 
pendence of the vibrational line width T{q) observed in 
many experiments on inelastic x-ray scattering in glasses. 
Therefore we think that it is necessary to take them into 
account in interpretation of experimental data. 
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